You are given the MODIS LAI data files for the year 2012 in the directory files/data for the UK (MODIS tile .
Read the LAI datasets into a masked array, using QA bit 0 to mask the data (i.e. good quality data only) and generate a movie of LAI.
You ought be getting familiar with this sort of problem, as it is very common.
Since we have the code in the main notes to read a single data dataset:
import gdal # Import GDAL library bindings
import numpy as np
# The file that we shall be using
# Needs to be on current directory
filename = 'files/data/MCD15A2.A2011185.h09v05.005.2011213154534.hdf'
g = gdal.Open(filename)
# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
print "Problem opening file %s!" % filename
else:
print "File %s opened fine" % filename
subdatasets = g.GetSubDatasets()
for fname, name in subdatasets:
print name
print "\t", fname
# Let's create a list with the selected layer names
selected_layers = [ "Lai_1km", "FparLai_QC" ]
# We will store the data in a dictionary
# Initialise an empty dictionary
data = {}
# for convenience, we will use string substitution to create a
# template for GDAL filenames, which we'll substitute on the fly:
file_template = 'HDF4_EOS:EOS_GRID:"%s":MOD_Grid_MOD15A2:%s'
for i, layer in enumerate ( selected_layers ):
this_file = file_template % ( filename, layer )
print "Opening Layer %d: %s" % (i+1, this_file )
g = gdal.Open ( this_file )
if g is None:
raise IOError
data[layer] = g.ReadAsArray()
print "\t>>> Read %s!" % layer
# scale the LAI
lai = data['Lai_1km'] * 0.1
# pull out the QC
qc = data['FparLai_QC']
# find bit 0
qc = qc & 1
# generate the masked array
laim = np.ma.array ( lai, mask=qc )
We can use this as a function that reads a single dataset:
def read_modis_lai(filename):
''' Read MODIS LAI data from filename
and return masked array
'''
g = gdal.Open(filename)
# g should now be a GDAL dataset, but if the file isn't found
# g will be none. Let's test this:
if g is None:
print "Problem opening file %s!" % filename
else:
print "File %s opened fine" % filename
subdatasets = g.GetSubDatasets()
#for fname, name in subdatasets:
#print name
#print "\t", fname
# Let's create a list with the selected layer names
selected_layers = [ "Lai_1km", "FparLai_QC" ]
# We will store the data in a dictionary
# Initialise an empty dictionary
data = {}
# for convenience, we will use string substitution to create a
# template for GDAL filenames, which we'll substitute on the fly:
file_template = 'HDF4_EOS:EOS_GRID:"%s":MOD_Grid_MOD15A2:%s'
for i, layer in enumerate ( selected_layers ):
this_file = file_template % ( filename, layer )
#print "Opening Layer %d: %s" % (i+1, this_file )
g = gdal.Open ( this_file )
if g is None:
raise IOError
data[layer] = g.ReadAsArray()
#print "\t>>> Read %s!" % layer
# scale the LAI
lai = data['Lai_1km'] * 0.1
# pull out the QC
qc = data['FparLai_QC']
# find bit 0
qc = qc & 1
# generate the masked array
laim = np.ma.array ( lai, mask=qc )
return laim
and then loop:
import glob
year = 2012
tile = 'h17v03'
files = np.sort(glob.glob('files/data/MCD15A2.A%d*.%s.*.hdf'%(year,tile)))
lai = []
for f in files:
lai.append(read_modis_lai(f))
# force it to be a masked array
lai = np.ma.array(lai)
# plot the data
import pylab as plt
# work out a consistent scaling
lai_max = np.max(lai)
for i,f in enumerate(files):
plt.figure(figsize=(7,7))
plt.imshow(lai[i],interpolation='none',vmin=0.,vmax=lai_max*0.75)
# remember filenames of the form
# files/data/MCD15A2.A2011185.h09v05.005.2011213154534.hdf'
file_id = f.split('/')[-1].split('.')[-5][1:]
print file_id
# plot a jpg
plt.title(file_id)
plt.colorbar()
plt.savefig('files/images/lai_uk_%s.jpg'%file_id)
# now make a movie ...
import os
cmd = 'convert -delay 100 -loop 0 files/images/lai_uk_*.jpg files/images/lai_uk2.gif'
os.system(cmd)